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GRAPHICAL  ABSTRACT 


•  Optical  images  show  that  GDL  fibers 
are  wavy  and  not  straight. 

•  Thermal  contact  resistance  (TCR)  de¬ 
creases  with  fiber  wavelength. 

•  TCR  increases  with  curvature,  diam¬ 
eter  and  amplitude  of  fiber  and  GDL 
porosity. 

•  TCR  does  not  change  with  fiber 
length. 

•  TCR  variations  with  geometric  pa¬ 
rameters  are  higher  at  lower 
compression. 
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A  new  analytical  model  is  developed  to  predict  the  thermal  contact  resistance  (TCR)  between  fibrous 
porous  media  such  as  gas  diffusion  layers  (GDLs)  of  polymer  electrolyte  membrane  fuel  cells  (PEMFCs)  and 
flat  surfaces  (bipolar  plates).  This  robust  model  accounts  for  the  salient  geometrical  parameters  of  GDLs, 
mechanical  deformation,  and  thermophysical  properties  of  the  contacting  bodies.  The  model  is  success¬ 
fully  validated  against  experimental  data,  and  is  used  to  perform  in  a  comprehensive  parametric  study  to 
investigate  the  effects  of  fiber  parameters  such  as  waviness  and  GDL  properties  on  the  TCR.  Fiber  waviness, 
diameter  and  surface  curvature,  as  well  as  GDL  porosity,  are  found  to  have  a  strong  influence  on  TCR 
whereas  fiber  length  does  not  affect  the  TCR  when  the  porosity  is  kept  constant.  Such  findings  provide 
useful  guidance  for  design  and  manufacturing  of  more  effective  GDLs  for  PEMFC  heat  management. 

The  analytic  model  can  be  readily  implemented  in  simulation  and  modeling  of  PEMFCs,  and  can  be 
extended  with  minor  modifications  to  other  fibrous  porous  media  such  as  fibrous  catalysts,  insulating 
media  and  sintered  metals. 

©  2014  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Durability,  reliability  and  stability  of  polymer  electrolyte 
membrane  fuel  cells  (PEMFCs)  are  strongly  dependent  on  heat  and 
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associated  water  management.  It  has  been  recently  shown  1,2] 
that  the  contact  resistance  between  the  PEMFC  components  may 
be  higher  than  their  bulk  resistances.  However,  contact  resistance 
has  typically  either  been  roughly  estimated  or  simply  overlooked 
in  performance  model  analyses  [3,4  .  This  is  mainly  due  to  the 
technical  barriers  and  challenges  in  measuring  and  predicting 
contact  resistance,  especially  over  a  range  of  compression,  see  e.g. 
Refs.  [5-12]. 


http://dx.doi.org/10.1016/jjpowsour.2014.04.149 
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Fig.  1.  (a)  GDL  surface  images  (Sigracet  GDLs)  and  (b)  Proposed  geometrical  modeling 
for  a  GDL 


Among  the  interfacial  resistances  in  PEMFCs,  the  contact  resis¬ 
tance  between  gas  diffusion  layer  (GDL)  and  neighboring  bipolar 
plate  (BPP)  is  of  special  interest  due  to  the  rib/channel  structure  of 
the  BPP  compressing  the  adjacent  GDL  The  complexity  of  the  GDL 
micro-structure  and  its  surface  morphology  have  lead  most  re¬ 
searchers  to  adopt  numerical  13,14]  and  experimental  [15-21] 
methods  rather  than  analytic  models.  In  our  previous  works 
(Sadeghifar  et  al.  [1,2,5  ),  the  thermal  contact  resistances  (TCR)  of 
different  GDLs  with  metallic  surfaces  and  also  with  graphite  BPP 
were  measured.  The  results  showed  that  the  contact  resistance  can 
be  as  high  as  the  GDL  bulk  resistance,  even  for  thick  GDLs  with 
thermal  conductivities  as  low  as  0.2  W  m_1  K_1.  A  few  numerical 
studies  have  also  been  performed  to  estimate  the  GDL-BPP  contact 
resistance.  Using  ANSYS  Fluent,  Nitta  et  al.  [16]  numerically  esti¬ 
mated  the  contact  resistance  between  a  GDL  and  a  graphite  plate, 
but  reported,  inconsistently  with  contact  mechanics  considerations 
and  available  experimental  observations,  that  the  GDL  thermal 
conductivity  is  independent  of  compression. 

The  objective  of  this  study  is  to  develop  a  mechanistic  analytic 
model  for  predicting  the  thermal  contact  resistance  between  a 


GDL-a  fibrous  porous  material-and  a  flat  surface,  e.g.  a  graphite  or 
metallic  BPP.  The  present  model  is  built  using:  i)  GDL  and  BPP 
salient  geometric  parameters,  such  as  waviness,  diameter,  distri¬ 
bution  and  orientation  of  fibers,  GDL  porosity;  ii)  applied  load, 
mechanical  deformation,  Hertzian  theory;  iii)  thermophysical 
properties  of  both  contacting  bodies,  i.e.,  fibrous  porous  medium 
and  flat  plate,  properties  such  as  thermal  conductivities  and 
effective  Young’s  modulus;  and  iv)  heat  conduction  in  GDL  fibres 
(spreading/constriction  resistances). 

2.  Model  development 

The  real  area  between  the  two  contacting  bodies  is  the  key 
parameter  in  determining  both  electrical  and  thermal  contact 
resistance.  The  contact  area  can  be  determined  using  geometrical 
and  mechanical  modeling.  The  TCR  at  the  interface  between  a 
fibrous  porous  medium  and  a  solid  flat  plate  can  be  obtained  using 
an  appropriate  thermal  model  that  includes  heat  transfer  in  both 
contacting  bodies  through  the  contacting  areas.  This  section  pro¬ 
vides  an  overview  of  different  elements  of  the  model. 

2.1.  Geometrical  model 

Based  on  images  of  different  GDLs,  the  fibers  are  assumed  to 
have  a  circular  cross-section,  as  was  the  case  in  our  previous  works 
[1,8  ,  see  Fig.  1.  In  almost  all  previous  GDL  geometric  models,  the 
fibers  are  assumed  to  be  straight,  i.e.  cylindrical.  When  modeling 
interfacial  phenomena,  which  are  highly  dependent  on  surface 
topography  and  morphology,  more  realistic  assumptions  are 
required.  Our  analysis  of  GDL  images  reveals  that  fibers  are  in  fact 
wavy,  as  can  be  seen  in  Fig.  2.  Consequently,  fibers  are  considered  as 
wavy  cylinders,  with  a  sinusoidal  profile  in  this  study,  see  Fig.  3.  The 
waviness  of  the  fibers  is  measured  optically  using  a  mechanical 
tester  microscope  (NANOVEA).  This  instrument  allows  micro-scale 
imaging/filming  while  scanning  the  sample  in  different  directions. 
The  statistical  data  for  fiber  waviness  and  amplitude  are  presented 
in  Fig.  4(a)  and  (b)  for  Sigracet  (SGL)  GDLs.  The  wavelength  (A)  and 
amplitude  (A)  data  fall  in  the  ranges  50-1900  pm  and  2-5df  for  the 
Sigracet  samples.  For  most  engineering  applications,  however,  us¬ 
ing  the  average  values  of  fiber  amplitude  and  wavelength  may  lead 
to  sufficient  accuracy. 


Fig.  2.  GDL  images  showing  waviness  of  the  fibers. 
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Fig.  3.  A  wavy  fiber  under  compression:  increasing  contact  areas  between  a  wavy  fiber  and  flat  surface  with  increasing  compression. 


The  proposed  model  is  developed  based  on  the  following  as¬ 
sumptions:  1)  steady  state  heat  transfer,  2)  constant  thermophys¬ 
ical  properties,  3)  wavy  GDL  fibers,  4)  smooth  flat  surface,  i.e.,  the 
surface  roughness  of  the  plate  is  neglected,  5)  no  out-of-flatness  on 
the  plate,  6)  static  mechanical  contact,  i.e.,  no  vibration  effects,  7) 
negligible  radiation  effects  due  to  low  typical  operating  tempera¬ 
ture  range  of  PEMFC  <100  °C  [1,2  ,  and  8)  first  loading  cycle  only, 
i.e.  no  hysteresis  effect  is  considered;  but  the  proposed  methodol¬ 
ogy  is  also  applicable  to  deformed  (cycled)  GDLs  if  the  deformed 
samples  geometric  parameters  are  available,  and  9)  contact  occurs 
in  a  vacuum  environment.  It  should  be  noted  that  the  effects  of  heat 
transfer  in  interstitial  gases  can  be  added  to  the  present  model 
using  the  model  of  Ref.  [6  . 

The  apparent  lengths  of  the  GDL  fibers  were  optically  measured 
and  their  lengths  were  obtained  by  calculating  the  length  of  the 
sinusoidal  arc.  Nevertheless,  because  of  the  small  waviness,  the  real 
length  can  be  simply  calculated  mathematically: 

If  =  2 (Ns  -  1)^A2  +A2/4  =  /fap  (1) 

The  total  number  of  fibers  in  one  layer  of  GDL,  Nft,  and  the 
number  of  contact  strips  that  each  fiber  can  form  on  the  flat  surface, 
Ns,  can  also  be  obtained  from  the  geometrical  data  of  the  GDL: 


on  the  contact  between  one  fiber  and  a  flat  surface  and,  then, 
extend  the  model  to  all  contacting  fibers. 


2.2.  Mechanical  model 


Thermal  energy  transfers  from  one  fiber  to  the  flat  surface 
through  the  contact  strips  at  the  interface.  The  resistance  to  heat 
conduction  depends  on  the  contact  area.  The  pressure  distribution 
can  be  calculated  for  an  elastic  flat  surface  in  contact  with  an  elastic 
wavy  cylinder,  as  shown  in  Fig.  3,  from  Ref.  [30  : 


where  x  is  a  space  variable  defined  along  the  cylinder  axis  and  p  is 
the  mean  pressure,  which  is  related  to  the  length  of  each  contact 
strip  using  Hertzian  theory  as: 


P 


sin' 


ira 

X 


(5a) 


Or  in  an  explicit  form  in  terms  of  contact  strip  length  [30]: 


4A(!  -e) 

TTdflf 


2  a 


21 

— arcsin 

7 r 


(5b) 


Ns=^+1  (3) 

The  geometrical  equations  and  GDLs  parameters  required  in  the 
model  are  presented  in  detail  in  Refs.  1,4,8]  and  are  simply  sum¬ 
marized  in  Table  1. 

Fig.  3  shows  the  ‘contact  strips’  in  the  contact  plane,  where  the 
contact  area  between  the  flat  surface  and  the  wavy  fibrous  porous 
medium  occurs.  The  GDL  fibers  have  a  relatively  small  amplitude  A 
and  a  wavelength  A.  The  dimensions  of  these  contact  strips,  a  and  t, 
grow  with  increasing  compression  as  illustrated  in  Fig.  3.  The 
contact  strips  dimensions  can  be  calculated  using  the  Hertzian 
contact  theory  and  by  applying  a  force  balance.  Here  we  first  focus 


Using  both  Hertzian  theory  and  force  balances  on  fibers  con¬ 
tacting  the  flat  surface,  the  width  of  each  contact  strip,  2 1,  can  be 
obtained  as  a  function  of  the  apparent  load  F  (the  load  applied  on 
the  entire  GDL,  see  Figs.  3  and  5): 

(  4Fdf  V 

2t  ~  yNf(Ns  -  1)aE*ir )  (6) 

It  is  evident  from  the  above  equations  that  the  geometrical 
parameters  of  the  contact  strips  are  highly  nonlinear  and  coupled, 
and  hence  no  explicit  functionality  can  be  derived  for  the  contact 
area  dimensions.  However,  for  most  applications,  the  value  of  a/A  is 
small,  and  the  term  sin(7ra/X)  reduces  to  7ra/X,  which  makes  the 
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Fig.  4.  Statistical  data  of  (a)  fiber  waviness  and  (b)  fiber  amplitude  for  Sigracet  (SGL) 
GDLs. 
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Having  calculated  the  contact  strip  length  (2a),  the  value  of  the 
strip  width  (2 1)  can  be  directly  obtained  from  Eq.  (6). 

GDL  surfaces  are  not  flat  and  show  a  random  distribution  of  the 
surface  asperities.  Following  Mikic  [24]  and  Bahrami  et  al.  [22  ,  we 
here  assume  a  Gaussian  distribution  for  the  distance  between  fibers 
and  the  flat  surface,  which  is  a  function  of  compression  as  shown 
schematically  in  Fig.  5.  Mikic  [24]  reported  a  relationship  for  the 
number  of  Gaussian  asperities  of  an  elastic  body  contacting  a  flat 
surface,  as  a  function  of  compression,  as: 


n  _  J_/m\2£xp(l:2Y^ 
s  16  v  a)  erfc(y) 


(8a) 


Y  =  erfc  1  (8b) 

where  m,  a,  y,  A ,  Hei,  P  are  asperity  slope,  GDL  effective  surface 
roughness,  dimensionless  mean  plane  separation,  apparent  (total) 
area,  elastic  micro-hardness,  and  pressure,  respectively,  see  also 
Table  1.  It  is  worth  mentioning  that  m  is  a  weak  function  of  a  (see 
Table  1)  and,  for  most  applications,  a  value  of  m  =  0.1  can  be  used 
[22].  The  number  of  fibers  contacting  the  surface  (Nf)  at  any 
compression  of  P  can  be  obtained  by  the  same  proportionality  as  Eq. 
(8)  proposes: 


Pi  =  jh 

Nft  nst 


where  N ft  is  the  total  number  of  fibers  as  given  in  able  1.  According 
to  Eq.  (9),  the  present  model  uses  only  the  proportionality  of  Eq.  (8), 
as  it  obtains  Nft  from  the  GDL  porosity.  Hence,  the  TCR  results  are 
not  sensitive  to  the  overall  GDL  roughness  with  a  range  of  l-3d/ 
reported  in  the  literature  [18,31,32]. 


analysis  simpler.  Hence,  for  most  applications  including  GDLs,  one 
can  obtain  the  length  of  each  contact  strip  (2a)  in  terms  of  acces¬ 
sible/measurable  parameters  (especially  F),  directly  from  the 
following  equation: 

n  _  1  {?)*{  2F 

Q  8tt4  ^ttA J  \Nf(Ns  -  \)E*df 


2.3.  Thermal  model 

Due  to  the  very  small  area  of  the  contact  strips,  the  heat 
transferred  from  one  fiber  to  the  flat  surface  encounters  a  large 
resistance,  mostly  referred  to  as  spreading/constriction  resistance 
[1,8  .  The  spreading/constriction  resistance  that  occurs  on  the  cyl¬ 
inder  (fiber)  side  can  be  obtained  by  Ref.  [33]: 


Table  1 

Summary  of  geometrical  and  mechanical  parameters  of  the  studied  GDLs  and  mechanical  variables  used  in  the  present  model. 


Symbol 

Parameter 

Units 

Value  or  equation 

Basis 

Eq./Fig. 

number 

E 

Young  modulus  of  fiber 
&  plate 

GPa 

210  &  210  [1,7,8] 

Meas. 

— 

u 

Poisson  ratio  of  fiber 
&  plate 

0.3  &  0.3  [1,7,8] 

Meas. 

k 

Thermal  conductivity  of 
fiber  &  plate 

W  m1  K1 

120  &  70  [1,4, 7, 8] 

Meas. 

— 

Ifap 

Apparent  fiber  length 

pm 

3000 

Meas. 

— 

df 

Fiber  diameter 

pm 

7.5  [1,8] 

Meas. 

— 

X 

Wavelength 

pm 

50-1900 

Meas. 

Fig.4(a) 

A 

Amplitude 

pm 

4  df 

Meas. 

Fig.4(b) 

A 

GDL  cross-sectional  area 
(apparent  surface  area) 

2 

m 

0.000507  [1,8] 

Meas. 

e 

Nominal  porosities  of 

SGL  24AA  &  25AA 

— 

0.88  &  0.92  [1,8] 

Meas. 

— 

a 

Roughness  of  GDLs 

SGL  24AA  &  25AA 

m 

17.0  ±  3.5  &  31.0  ±  4.5 

[18,1] 

Meas. 

— 

m 

Asperities  slope  for  GDL 

— 

0.076g°’52  [22-29] 

Calc. 

— 

He  1 

Effective  elastic  modulus 

Pa 

Ej?  [22-24] 

Calc. 

— 
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Fig.  5.  (a)  Gaussian  distribution  of  the  fibers  at  the  GDL  surface  and  (b)  the  variations  in  the  number  and  size  of  contact  areas  with  respect  to  compression  (F"  >  F  >  F). 
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And  that  on  the  flat  surface  is  given  as  [33]: 
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where  u  is  the  half  width  of  the  rectangular  heat  channels  on  which 
the  isoflux  (rectangular)  contact  strips  are  located  centrally,  which 
can  be  obtained  by: 


2  u 


A 

XNf(Ns  -  1) 


(12) 


It  is  worth  mentioning  that  the  summation  of  all  the  heat  chan¬ 
nels  on  the  flat  surface  should  be  equal  to  the  GDL  cross-sectional 
area.  A  detailed  explanations  on  the  heat  channel  concept  for  con¬ 
tact  resistance  estimation  can  be  found  elsewhere  !  22,33]. 

The  TCR  for  each  contact  strip  is  the  sum  of  the  two  resistances 
given  by  Eqs.  (10)  and  (11).  Ultimately,  the  TCR  between  the  GDL 
and  the  flat  surface  can  be  obtained  using  a  parallel  summation  of 
all  the  contact  strip  resistances: 


Rc  + 

Nf(Ns-  1) 


(ii) 


TCR 


(13) 
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The  TCR  averaged  over  all  the  fiber  waviness  can  be  obtained  by 
the  summation  of  the  TCRs  for  different  waviness  based  on  their 
occurrence  probability  for  the  GDL  of  interest: 

TCRaVeA  =  EPAjTCRAj-  (14) 

j=  1 

where  pjj  is  the  probability  of  having  the  waviness  of  Xj,  which  is 
given  in  Fig.  4a,  and  TCR^j  is  the  value  of  TCR  obtained  from  Eq.  (13) 
for  fiber  waviness  of  Xj.  It  should  be  noted  that  since  the  majority  of 
the  fiber  amplitudes  usually  lies  within  a  narrow  range,  using  an 
average  value  of  the  amplitude  is  appropriate  for  the  TCR 
calculations. 

3.  Results  and  discussion 

Once  the  geometrical  parameters  of  a  GDL  are  determined,  one 
can  calculate  the  associated  TCR  as  a  function  of  compression.  A 
computer  code  is  developed  in  MATLAB  to  facilitate  the  TCR  =  f(P,  A, 
A ,  e,  df,  dfsc,  k)  calculations  for  the  parametric  study.  It  should  be 
noted  that  fiber  specifications  (A,  A ,  df,  dfso  If)  and  porosity  (e)  of  the 
virgin  GDLs  are  geometric  parameters  required  as  input  for  the 
model  and  are  hence  used  in  the  parametric  study.  The  effect  of 
these  parameters  on  TCR  is  accounted  for  in  conjunction  with 
compression.  For  instance,  with  increasing  compression,  the 
number  of  fibers  contacting  the  plate  and  the  contact  area  di¬ 
mensions  (a  and  t)  increase,  as  explained  in  Section  “Mechanical 
model”  and  shown  schematically  in  Figs.  3  and  5.  The  wavelength 
however  does  not  change.  It  should  be  noted  that  all  the  TCR  cal¬ 
culations  in  this  study  are  based  on  one-inch  circular  GDLs 
(A  =  0.000507  m2).  The  TCR  per  unit  area  of  a  GDL  (TCRa)  can  be 
obtained  as  TCRa  =  TCR  x  A. 

3.1.  Model  validation 

Fig.  6  compares  the  results  of  the  present  model  with  our 
experimental  data  [1  for  two  GDLs  SGL  24AA  and  SGL  25AA.  The 
results,  shown  by  solid  lines,  are  in  good  agreement  with  the 
experimental  data  and  the  model  captures  the  experimental  TCRs 
over  a  wide  range  of  compression.  Fig.  6  shows  that  the  TCR  de¬ 
creases  with  increasing  compression,  but  the  slope  becomes  shal¬ 
lower  and  eventually,  at  high  compression,  the  impact  of 
compression  on  TCR  becomes  negligible,  as  expected.  For  SGL  GDLs 
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Fig.  7.  Impact  of  porosity  on  the  TCR  at  different  compressions. 

having  an  average  amplitude  of  A/df  =  4  (Fig.  4b),  the  model  is  in 
good  agreement  with  experiments,  except  for  very  low  pressures 
(2  bar)  where  the  model  provided  less  accurate  estimations  of  the 
TCR  at  the  complex  interface  of  the  GDL-plate. 

3.2.  Parametric  study 

In  order  to  investigate  the  effects  of  key  parameters  on  the  GDL- 
BPP  TCR,  a  parametric  study  is  performed  in  this  section.  For  each 
case  study,  when  a  parameter  is  changed,  the  other  parameters  are 
kept  constant  unless  otherwise  mentioned. 

3.2.1.  GDL  porosity  (e) 

The  most  important  and  experimentally  accessible  parameter  of 
a  porous  medium  is  porosity,  whose  effect  on  all  the  transport 
properties  is  usually  significant.  The  effect  of  porosity  on  the  TCR  is 
shown  in  Fig.  7  for  three  compressions  of  4,  10  and  20  bar.  The 
curves  show  similar  trends  for  all  three  compressions  and  with 
increasing  porosity,  the  TCR  increases  as  well.  One  important  point 
is  that  for  each  curve,  especially  at  lower  compressions,  there  is  a 
specific  value  of  the  porosity  at  which  a  noticeable  increase  can  be 
observed  at  the  rate  of  the  TCR  variations.  This  specific  porosity, 
determined  here  at  the  average  slope  of  each  TCR-s  curve  on  Fig.  7, 
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Fig.  6.  Comparison  of  model  results  (solid  lines)  and  experimental  data  for  SGL24AAand  25AA  (the  GDL  specifications  has  already  been  given  in  Fig.  4  and  Table  1:  df  =  7.5  pm  [1,8], 
If  =  3000  pm,  A/df  =  4,  A  =  0.000507  m2  [1,8],  e  =  88.3  &  92.1%  for  the  GDLs  SGL  24AA  &  25AA,  respectively). 
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is  approximately  89%  at  all  three  compressions.  This  point  should 
be  accounted  for  in  GDL  manufacturing  and  fuel  cell  design. 

3.2.2.  Fiber  length  (If) 

Fig.  8  shows  that  the  fiber  length  does  not  have  any  effect  on  the 
TCR  at  different  constant  porosities.  This  is  due  to  a  trade-off  be¬ 
tween  the  number  of  fibers  (Nf)  and  the  number  of  contact  strips 
each  fiber  can  form  on  the  surface  (Ns)  at  a  constant  porosity. 
Increasing  fiber  lengths  at  constant  porosity  is  equivalent  to 
decreasing  the  number  of  fibers;  however,  the  number  of  contact 
strips  increases  as  well,  so  that  the  total  number  of  strips  remains 
constant  throughout  the  interface.  As  a  result,  the  contact  area  does 
not  change  with  changing  average  fiber  length  at  a  constant 
porosity,  as  observed  in  Fig.  8. 

3.2.3.  Fiber  wavelength  (X) 

One  of  the  two  key  parameters  in  determining  the  GDL-BPP  TCR 
is  the  fiber  wavelength  and  amplitude,  which  determines  the  de¬ 
gree  of  waviness  that  the  GDL  fibers  have.  The  impact  of  fiber 
wavelength  is  presented  in  Fig.  9  for  high  and  low  porosity.  TCR 
decreases  noticeably  with  increasing  fiber  wavelength  with  a  slope 
that  is  much  more  pronounced  for  lower  compression  (P  =  4  bar) 
and  shorter  wavelengths.  The  rate  at  which  the  TCR  decreases 
become  markedly  lower  where  A  reaches  a  specific  wavelength, 
determined  at  the  average  slope  of  any  of  the  TCR  curves  in  Fig.  9, 
and  corresponds  approximately  to  600  and  700  pm  for  the  low  and 
high  porosity  cases,  respectively.  This  point  can  be  an  important 
consideration  from  the  viewpoint  of  GDL  manufacturing. 

3.2.4.  Fiber  amplitude  (A  ) 

Fig.  10  shows  that  the  fibre  amplitude  can  have  a  strong  effect  on 
TCR.  However,  for  typical  fiber  amplitudes,  i.e.,  3-4df,  the  effect  on 
TCR  is  not  significant.  At  larger  amplitudes,  especially  higher  than 
3df,  this  effect  becomes  negligible  at  medium  to  high  compression. 
A  similar  trend  was  already  seen  in  Fig.  9  regarding  the  effect  of  the 
fiber  wavelength  on  TCR. 

Comparing  the  results  of  zero  amplitude  (A  =  0  pm),  which 
corresponds  to  straight  fibers,  with  those  of  other  amplitudes 
clearly  shows  that  the  fibers  waviness  can  increase  the  TCR  dras¬ 
tically,  even  though  its  effect  on  the  TCR  drops  noticeably  beyond  a 
specific  value  of  approximately  1.5df.  This  point  should  be  taken 
into  account  in  design  of  both  GDLs  and  fuel  cell  stacks  in  terms  of 
heat  management.  For  instance,  GDLs  consisting  of  very  wavy  fibers 
are  expected  to  have  a  large  contact  resistance  with  the 


Fig.  8.  Effect  of  fiber  length  on  the  TCR  at  three  compressions  of  4, 10  and  20  bar. 


Fig.  9.  Impact  of  fiber  wavelength  on  the  TCR  for  two  low  and  high  porosities  at  three 
different  compressions. 


neighboring  BPP  and,  hence,  may  be  less  effective  in  PEMFCs  in 
terms  of  heat  (and  electrical)  management. 

3.2.5.  Fiber  diameter  (df) 

Fig.  11  shows  that  fiber  diameter  can  significantly  change  the 
GDL-BPP  TCR  at  low  contact  pressures  and  with  increasing 
compression,  its  effect  decreases.  At  high  compressions,  for  a 
typical  range  of  GDL  fiber  diameter,  i.e.,  7-10  pm  [1,8],  the  effect  of 
fiber  diameter  on  TCR  becomes  insignificant. 

It  is  worth  noting  that  the  functionality  of  the  TCR  with  fiber 
diameter  is  practically  linear  and,  hence,  no  critical  value  for  fiber 
diameter  can  be  determined  with  respect  to  the  TCR.  However,  as 
Fig.  11  shows,  with  smaller  fiber  diameters,  the  TCR  decreases  and, 
eventually,  converges  to  a  very  low  value  of  0.025  K  W-1,  inde¬ 
pendent  of  the  compression.  This  strongly  suggests  that  GDL 
manufacturers  explore  production  of  GDLs  with  thinner  fibers  as 
long  as  other  factors  such  as  manufacturing  capabilities  or  GDL 
mechanical  strength  allow. 

3.2.6.  Fiber  surface  curvature  at  the  contact  interface  (dfSC) 

Some  types  of  fibers  are  not  completely  cylindrical;  e.g.,  they 
may  be  ellipsoidal,  while  having  the  same  volume  as  a  cylindrical 
fiber.  The  deviation  from  a  perfectly  cylindrical  shape  can  be 
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Fig.  11.  Effect  of  fiber  diameter  on  TCR  at  three  compressions  of  4, 10,  and  14  bar:  Due 
to  almost  linear  functionality  of  TCR  with  fiber  diameter,  no  critical  value  for  fiber 
diameter  can  be  determined.  However,  the  TCR  reduces  to  0.025  K  W  1  as  fiber 
diameter  approaches  zero  (all  the  other  parameters  except  for  the  number  of  fibers  are 
kept  constant  to  keep  the  GDL  porosity  unchanged  (e  =  0.88)). 

accounted  for  by  considering  the  effect  of  fiber  surface  curvature  at 
the  contact  interface.  The  effect  on  TCR  of  fiber  curvature  (which  in 
the  limit  of  a  cylinder  is  the  fiber  diameter)  can  be  analyzed  without 
changing  any  parameters  to  keep  the  GDL  porosity  constant  (e.g., 
without  changing  the  number  of  fibers).  Fig.  12  shows  that,  whereas 
TCR  varies  linearly  with  fiber  diameter,  its  dependence  on  fiber 
surface  curvature  is  highly  non-linear  and  is  particularly  significant 
at  very  small  curvatures.  Comparing  Figs.  11  and  12  also  shows  that, 
for  typical  values  of  fiber  diameter  and  fiber  curvature,  the  TCR  is 
not  as  sensitive  to  the  fiber  surface  curvature  as  to  the  fiber 
diameter  (which  determines  both  fiber  surface  curvature  and  fiber 
volume).  Overall,  the  points  presented  in  this  study  can  also  be 
taken  into  account  in  determining  the  TCR  of  GDLs  with  catalyst 
coated  membranes  (CCMs)  [34]  and  microporous  layers  (MPLs) 
[35]. 

4.  Summary  and  conclusion 

A  new  and  robust  mechanistic  model  was  developed  for  pre¬ 
dicting  the  thermal  contact  resistance  between  a  general  fibrous 
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Fig.  12.  Effect  of  fiber  surface  curvature  at  the  contact  interface  on  the  TCR:  df  =  dfSC  in 
all  the  equations  except  for  Eq.  (2)  which  relates  the  GDL  porosity  with  the  number 
and  volume  of  the  fibers. 


medium  and  a  flat  surface,  and  was  applied  to  analyze  the  contact 
between  PEMFC  gas  diffusion  layers  and  bipolar  plates.  The  model 
accounts  for  the  salient  geometrical  parameters,  mechanical 
deformation  and  thermophysical  properties,  and  captures  accu¬ 
rately  experimentally-observed  behavior  over  a  wide  range  of 
compression,  except  for  very  low  compression  (less  than  ~2  bar) 
where  the  model  provided  rough  predictions  of  the  TCR. 

The  model  allows  systematic  investigation  of  the  impact  of  the 
main  GDL  micro  structural  properties  on  TCR  without  having  to  rely 
entirely  on  more  difficult  experimental  measurements.  The  simple 
and  general  model  developed  in  this  study  can  be  readily  imple¬ 
mented  into  fuel  cell  models  and  may  also  be  used  for  modeling  the 
behavior  of  any  other  fibrous  porous  materials  such  as  fibrous 
insulations. 
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Nomenclature 

A:  GDL  cross-sectional  area,  m2 

a:  half  of  contact  strip  length,  m 

b:  length  of  non-contacting  area,  b  =  7/2 -a,  m 

Calc.:  based  on  calculation  or  derivation 

df  fiber  diameter,  m 

dfSC:  diameter  of  fiber  surface  curvature  at  the  contact  interface,  m 
E:  Young  (elastic)  modulus,  Pa 

E*:  =  (1  -  y2/Ec  +  1  -  vj/EF)  Effective  Young  (elastic)  modulus,  Pa 

Exp.:  experimental  value 

F:  force  applied  on  the  entire  GDL,  N 

GDL:  gas  diffusion  layer 


k:  thermal  conductivity,  W  m-1  K-1 
If:  fiber  length,  m 
lfap:  apparent  fiber  length,  m 
Meas.:  based  on  measurement 

Nf.  number  of  fibers  contacting  the  surface  at  the  pressure  of  P 

Nft:  total  number  of  fibers  contacting  the  surface 

Ns:  number  of  troughs  (strips)  each  fiber  has  (forms  on  the  flat  surface) 

ns:  number  of  asperities  contacting  the  surface  of  another  body  in  contact 

nst:  total  number  of  asperities  contacting  the  surface 

P:  pressure  on  the  entire  GDL,  Pa 

p:  mean  pressure,  Pa 

PEMFC:  polymer  electrolyte  membrane  fuel  cell 
R:  thermal  resistance,  K  W-1 

Rc:  spreading/constriction  resistance  on  cylinder  side  resistance,  K  W-1 

RF:  spreading/constriction  resistance  on  flat  surface  side,  K  W-1 

t:  half  of  strip  width,  m 

TCR:  thermal  contact  resistance,  K  W-1 

TCRa:  specific-area  TCR:  Thermal  contact  resistance  per  unit  area  (=TCR  x  A), 
K  W1  m2 

u:  half  of  the  width  of  the  rectangular  channel  area,  m 
x:  a  space  variable  defined  along  the  cylinder  axis 

Greek  letter 

A:  amplitude  of  fiber  waviness,  m 
y:  non-dimentional  separation 
e:  porosity 

A:  wavelength  of  fiber  waviness,  m 

a:  root  mean  square  (RMS)  roughness  of  GDL  (=df  for  the  parametric  study),  m 
v:  Poisson’s  ratio 

Subscript 

C:  cylinder 
f:  fiber 
F:  flat  surface 
m:  measured 
n;  summation  index 
t:  total  value 


